/**
 * 预处理文件: data.txt  第1: int nSize(矩阵中非零个数)
 *N个int值表示每行中第一个数值的偏移 nSize个列 nSize个double类型的值
 *v1_b.txt 下面N个double类型的值
 *v1_x0.txt下面N个double类型的值
 *
**/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>
#include <omp.h>

#define Line 19
typedef struct
{
  int nx;
  int ny;
  int nz;
  int myrank;
}Info;

double *AV,*MV,*b,*x0,*diag;
int  *AJ;
int nSize,N,nN;  //nSize=NX*NY*NZ   N=nx*ny*nz nN=总长度 

Info InitInfo(int NX,int NY,int NZ,int PX,int PY,int PZ)
{
   Info result;
   int nx,ny,nz;
   nx=NX/PX;  ny=NY/PY; nz=NZ/PZ; 
   result.nx=NX/PX;  result.ny=NY/PY;  result.nz=NZ/PZ; 
   int myrank;
   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
  
   result.myrank=myrank;
   N=nx*ny*nz; 
   nN=N+2*nx*nz+2*ny*nz+4*nz;  

   return result;
}

void LoadData(Info *info,char *filename_A,char *filename_x0,char *filename_b)
{
  //printf("****-->myrank,nx,ny,nz  %d %d %d %d\n",info->myrank,info->nx,info->ny,info->nz);
    int nx,ny,nz,n;
    int ind=0;
    int i,j,k,m_j;
    int myrank=info->myrank;
    int s_ab,s_bc,s_cd ,s_ad,s_a,s_b,s_c,s_d;
    nx=info->nx;  ny=info->ny;  nz=info->nz;
 
    s_ab=N;   s_bc=s_ab+ny*nz; 
    s_cd=s_bc+nx*nz;
    s_ad=s_cd+ny*nz;
    s_a=s_ad+nx*nz;
    s_b=s_a+nz;
    s_c=s_b+nz;
    s_d=s_c+nz; 
    
    AV=(double*)malloc(N*Line*sizeof(double));
    AJ=(int*)malloc(N*Line*sizeof(int));
    x0=(double*)malloc(nN*sizeof(double));
    b=(double*)malloc(N*sizeof(double));  
    
    long bSeek=(myrank/2)*2*N;
    int bstep=ny*nz;
    FILE *fp2=fopen(filename_b,"rb");
    if(fp2==NULL){
        printf("can't open read file\n");
        exit(0);
    }
    fseek(fp2,bSeek*8,SEEK_CUR);
    ind=0;
    if(myrank%2)
    {
      for(i=0;i<nx;i++){
        ind=i*bstep;
        fseek(fp2,bstep*8,SEEK_CUR);
        fread(&b[ind],8,bstep,fp2);
      }
    }else{
      for(i=0;i<nx;i++){
        ind=i*bstep;
        fread(&b[ind],8,bstep,fp2);
        fseek(fp2,bstep*8,SEEK_CUR);
      }
   }
   fclose(fp2);

    FILE *fp3=fopen(filename_x0,"rb");
    if(fp3==NULL){
        printf("can't open read file\n");
        exit(0);
    }
    double *mx0=(double *)malloc(nSize* sizeof(double));
    fread(mx0,8,nSize,fp3);
    fclose(fp3);
//    printf("------>%.30lf\n",b[1]); 
   
    FILE *fp1=fopen(filename_A,"rb");
    if(fp1==NULL){
        printf("can't open read file\n");
        exit(0);
    }

    int size;
    fread(&size,4,1,fp1);
    int *mAJ=(int*)malloc(size*sizeof(int));
    fread(mAJ,4,size,fp1);
    
    long seekSize=(myrank/2)*2*N*Line;  //需要step长度
    long stepSize=ny*nz*Line;
    int aj_ind=seekSize;   //AJ对应位置     
    fseek(fp1,seekSize*8,SEEK_CUR);
 
    ind=0;
    if(myrank%2){    //基数块
      aj_ind=aj_ind+stepSize;
      for(i=0;i<nx;i++){
  ind=stepSize*i;  //storage location
        fseek(fp1,stepSize*8,SEEK_CUR);
        fread(&AV[ind],8,stepSize,fp1); 
     }
    }else{
      for(i=0;i<nx;i++){
  ind=stepSize*i;  //storage location
        fread(&AV[ind],8,stepSize,fp1); 
        fseek(fp1,stepSize*8,SEEK_CUR);
     }
   }
   fclose(fp1);
   //printf("----aj_ind--->%d\n",aj_ind);

   ind=0;
   int m_ind=0;
   int b_ind=0;
    for(i=0;i<nx;i++)
    {
        for(j=0;j<ny;j++){
            for(k=0;k<nz;k++){
                for(m_j=1;m_j<=Line;m_j++){
                switch (m_j)
                {
                    case 1:
                    {
                        m_ind=i*ny*nz+j*nz+k;
                        AJ[ind]=m_ind; x0[m_ind]=mx0[mAJ[aj_ind]];
                        aj_ind++;   ind++;
                        break;
                    }
                    case 2:
                    {
                        m_ind=(i-1)*ny*nz+j*nz+k;
                        if(i==0){
                            m_ind=s_ab+j*nz+k; 
                        }
                        AJ[ind]=m_ind; x0[m_ind]=mx0[mAJ[aj_ind]];
                        aj_ind++;   ind++;
                        break;
                    }
                    case 3:
                    {
                        m_ind=(i+1)*ny*nz+j*nz+k;
                        if(i==(nx-1)){
                            m_ind=s_cd+j*nz+k;
                        }
                        AJ[ind]=m_ind; x0[m_ind]=mx0[mAJ[aj_ind]];
                        aj_ind++;   ind++;
                        break;
                    }
                    case 4:
                    {
                        m_ind=i*ny*nz+(j-1)*nz+k;
                        if(j==0){
                            m_ind=s_ad+i*nz+k;
                        }
                        AJ[ind]=m_ind; x0[m_ind]=mx0[mAJ[aj_ind]];
                        aj_ind++;   ind++;
                        break;
                    }
                    case 5:
                    {
                        m_ind=i*ny*nz+(j+1)*nz+k;
                        if(j==(ny-1)){
                            m_ind=s_bc+i*nz+k; 
                        }
                        AJ[ind]=m_ind; x0[m_ind]=mx0[mAJ[aj_ind]];
                        aj_ind++;   ind++;
                        break;
                    }
                    case 6:
                    {
                        m_ind=(i+1)*ny*nz+(j+1)*nz+k;
                        if(i==(nx-1)){
                            m_ind=s_cd+(j+1)*nz+k; 
                        }
                        if(j==(ny-1)){
                            m_ind=s_bc+(i+1)*nz+k;
                        }
                        if((i==nx-1) && (j==(ny-1)))
                        {   
                            m_ind=s_c+k;
                        }
                        AJ[ind]=m_ind; x0[m_ind]=mx0[mAJ[aj_ind]];
                        aj_ind++;   ind++;
                        break;
                    }
                    case 7:
                    {
                        m_ind=(i+1)*ny*nz+(j-1)*nz+k;
                        if(i==(nx-1)){
                            m_ind=s_cd+(j-1)*nz+k; 
                        }
                        if(j==0){
                            m_ind=s_ad+(i+1)*nz+k;
                        }
                        if((i==nx-1) && (j==0))
                        {   
                            m_ind=s_d+k;
                        }
                        AJ[ind]=m_ind; x0[m_ind]=mx0[mAJ[aj_ind]];
                        aj_ind++;   ind++;
                        break;
                    }
                    case 8:
                    {
                        m_ind=(i-1)*ny*nz+(j-1)*nz+k;
                        if(i==0){
                            m_ind=s_ab+(j-1)*nz+k; 
                        }
                        if(j==0){
                            m_ind=s_ad+(i-1)*nz+k;
                        }
                        if((i==0) && (j==0))
                        {   
                            m_ind=s_a+k;
                        }
                        AJ[ind]=m_ind; x0[m_ind]=mx0[mAJ[aj_ind]];
                        aj_ind++;   ind++;
                        break;
                    }
                    case 9:
                    {
                        m_ind=(i-1)*ny*nz+(j+1)*nz+k;
                        if(i==0){
                            m_ind=s_ab+(j+1)*nz+k; 
                        }
                        if(j==(ny-1)){
                            m_ind=s_bc+(i-1)*nz+k;
                        }
                        if((i==0) && (j==(ny-1)))
                        {   
                            m_ind=s_b+k;
                        }
                        AJ[ind]=m_ind; x0[m_ind]=mx0[mAJ[aj_ind]];
                        aj_ind++;   ind++;
                        break;
                    }
                    case 10:
                    {
                        m_ind=i*ny*nz+j*nz+k-1;
                        if(k==0){
                            m_ind=N+1; 
                        }else{
         x0[m_ind]=mx0[mAJ[aj_ind]];
      }
                        AJ[ind]=m_ind; 
                        aj_ind++;   ind++;
                        break;
                    }
                    case 11:
                    {
                        m_ind=(i-1)*ny*nz+j*nz+k-1;
                        if(k==0){
                            m_ind=N+1;
                        }
                        else{
                            if(i==0){
                                m_ind=s_ab+j*nz+k-1;
                            }
          x0[m_ind]=mx0[mAJ[aj_ind]]; 
                        }
                        AJ[ind]=m_ind;
                        aj_ind++;   ind++;
                        break;
                    }
                    case 12:
                    {
                        m_ind=(i+1)*ny*nz+j*nz+k-1;
                        if(k==0){
                            m_ind=N+1;
                        }
                        else{
                            if(i==(nx-1)){
                                m_ind=s_cd+j*nz+k-1;
                            }
      x0[m_ind]=mx0[mAJ[aj_ind]]; 
                        }
                        AJ[ind]=m_ind; 
                        aj_ind++;   ind++;
                        break;
                    }
                    case 13:
                    {
                        m_ind=i*ny*nz+(j-1)*nz+k-1;
                        if(k==0){
                            m_ind=N+1;
                        }
                        else{
                            if(j==0){
                                m_ind=s_ad+i*nz+k-1;
                            }
       x0[m_ind]=mx0[mAJ[aj_ind]]; 
                        }
                        AJ[ind]=m_ind; 
                        aj_ind++;   ind++;
                        break;
                    }
                    case 14:
                    {
                        m_ind=i*ny*nz+(j+1)*nz+k-1;
                        if(k==0){
                            m_ind=N+1;
                        }
                        else{
                            if(j==(ny-1)){
                                m_ind=s_bc+i*nz+k-1;
                            }
         x0[m_ind]=mx0[mAJ[aj_ind]]; 
                        }
                        AJ[ind]=m_ind; 
                        aj_ind++;   ind++;
                        break;
                    }
                    case 15:
                    {
                        m_ind=i*ny*nz+j*nz+k+1;
                        if(k==(nz-1)){
                            m_ind=N+1;
                        }else{
       x0[m_ind]=mx0[mAJ[aj_ind]];
      }
                        AJ[ind]=m_ind; 
                        aj_ind++;   ind++;
                        break;
                    }
                    case 16:
                    {
                        m_ind=(i-1)*ny*nz+j*nz+k+1;
                        if(k==(nz-1)){
                            m_ind=N+1;
                        }
                        else{
                            if(i==0){
                                m_ind=s_ab+j*nz+k+1;
                            } 
           x0[m_ind]=mx0[mAJ[aj_ind]];
                        }
                        AJ[ind]=m_ind;
                        aj_ind++;   ind++;
                        break;
                    }
                    case 17:
                    {
                        m_ind=(i+1)*ny*nz+j*nz+k+1;
                        if(k==(nz-1)){
                            m_ind=N+1;
                        }
                        else{
                            if(i==(nx-1)){
                                m_ind=s_cd+j*nz+k+1;
                            }
      x0[m_ind]=mx0[mAJ[aj_ind]]; 
                        }
                        AJ[ind]=m_ind;
                        aj_ind++;   ind++;
                        break;
                    }
                    case 18:
                    {
                        m_ind=i*ny*nz+(j-1)*nz+k+1;
                        if(k==(nz-1)){
                            m_ind=N+1;
                        }
                        else{
                            if(j==0){
                                m_ind=s_ad+i*nz+k+1;
                            }
          x0[m_ind]=mx0[mAJ[aj_ind]]; 
                        }
                        AJ[ind]=m_ind; 
                        aj_ind++;   ind++;
                        break;
                    }
                    case 19:
                    {
                        m_ind=i*ny*nz+(j+1)*nz+k+1;
                        if(k==(nz-1)){
                            m_ind=N+1;
                        }
                        else{
                            if(j==(ny-1)){
                                m_ind=s_bc+i*nz+k+1;
                            }
         x0[m_ind]=mx0[mAJ[aj_ind]]; 
                        }
                        AJ[ind]=m_ind;
                        aj_ind++;   ind++;
                        break;
                    }
                    default:
                        break;
                }
        }
            }
          }
      aj_ind=aj_ind+stepSize;
    }
   free(mAJ); free(mx0); 
}

void getAmVector(double *result,double *data) //矩阵向量相乘
{
    int ind=0;
    int i,j;
    double tmp;
    #pragma omp parallel for private(i,j,tmp,ind) shared(AJ,data,AV,data,result)
    for(i=0;i<N;i++) //N-1
    {    
        tmp=0.0;
        ind=i*Line;
        for(j=0;j<Line;j++){
           tmp=tmp+AV[ind+j]*data[AJ[ind+j]];
        }
        result[i]=tmp;
    }
}

double VectorDianJi(double *dataA,double *dataB)
{
    int i;
    double sum=0.0;
   #pragma omp parallel for private(i) shared(dataA,dataB) reduction(+:sum)
    for(i=0;i<N;i++){
        sum=sum+dataA[i]*dataB[i];
    }
    return sum;
}


void JacobiCal(double *Rm,double *R) //矩阵Rm=M^(-1)*R Jacobi迭代
{
    int ind=0;
    int i,j;
    double tmp,m_rst;
    #pragma omp parallel for private(i,j,tmp,ind,m_rst) shared(AV,AJ,R,Rm)
    for(i=0;i<N;i++){
        ind=i*Line;
        tmp=0.0;    m_rst=0.0;
        for(j=0;j<Line;j++){
            if(AJ[ind+j]==i){  //对角线上的元素
                tmp=AV[ind+j];
            }else{
    m_rst=m_rst+AV[ind+j]*R[AJ[ind+j]];
            }
       }
       Rm[i]=(R[i]-m_rst)/tmp;
    }
}

double getjk(int j,int k)//获取Ajk值
{
   int i,ind;
   double reVal=0.0;
   ind=j*Line;
   for(i=0;i<Line;i++)
   {
      if(AJ[ind+i]==k){
         reVal=MV[ind+i];
         break;     
       }      
   }
   return reVal;
}

void getLUmatrixA() 
{
    diag=(double*)malloc(N*sizeof(double));   
    MV=(double *)malloc(N*Line* sizeof(double));
    memcpy(MV,AV,N*Line*sizeof(double));
    int i,j,k,ind;
    for(i=0;i<N;i++)
    {  
      ind=i*Line;
      for(j=0;j<Line;j++){
          if(AJ[ind+j]==i){
             diag[i]=MV[ind+j];  //对角元素
             break;
          }
      }
    }
    ind=0;
    for(i=0;i<N;i++){
      ind=i*Line; 
      for(j=0;j<Line;j++)
      { 
         if(AJ[ind+j]<i){
            MV[ind+j]=MV[ind+j]/diag[AJ[ind+j]];
            for(k=0;k<Line;k++){
               if((AJ[k+ind]>AJ[ind+j]) && (AJ[k+ind]<N))
               {
                  MV[ind+k]=MV[ind+k]-MV[ind+j]*getjk(AJ[ind+j],AJ[ind+k]);
                }
            }
        }
      }
  }
}


void getMatrixLU(double *Rm,double *R) 
{
   int i,j,ind;
   double *z=(double*)malloc(N*sizeof(double)); 
   double tmp;

   //#pragma omp parallel for private(i,j,tmp,ind) shared(MV,AJ,R,z)
   for(i=0;i<N;i++)
   {
     ind=i*Line;
     for(j=0;j<Line;j++){
          if(AJ[ind+j]==i){
             diag[i]=MV[ind+j];  
             break;
          }
      }
       tmp=0.0;
       for(j=0;j<Line;j++){
         if(AJ[ind+j]<i){
            tmp=tmp+MV[ind+j]*z[AJ[ind+j]];
         }
      }
      z[i]=R[i]-tmp/diag[i];
   }

 //  #pragma omp parallel for private(i,j,tmp,ind) shared(MV,AJ,Rm,z,diag)
   for(i=N-1;i>=0;i--)
   {
       ind=i*Line;
       tmp=0.0;
       for(j=0;j<Line;j++){
         if((AJ[ind+j]>i) && (AJ[ind+j]<N)){
            tmp=tmp+MV[ind+j]*Rm[AJ[ind+j]];
         }
       }
      Rm[i]=(z[i]-tmp)/diag[i];
   }
     free(z);
}


//get send data
void GetSendData(double *sendAB,double *sendBC,double *sendCD,double *sendAD,double *sendA,double *sendB,double *sendC,double *sendD,double *data,Info *info)
{
    int nx,ny,nz;
    nx=info->nx;        ny=info->ny;    nz=info->nz;
    
    int i,j,k;
    int ind=0,x_ind;
   #pragma omp parallel for private(j,k,ind,x_ind) shared(data,ny,nz) 
    for(j=0;j<ny;j++){
      for(k=0;k<nz;k++){
          ind=j*nz+k;
          x_ind=j*nz+k;
          sendAB[ind]=data[x_ind];
          x_ind=(nx-1)*ny*nz+j*nz+k;
          sendCD[ind]=data[x_ind];
      }
    }   
   // ind=0;
    x_ind=0;
   #pragma omp parallel for private(j,k,ind,x_ind) shared(data,ny,nz) 
    for(i=0;i<nx;i++){
      for(k=0;k<nz;k++){
   ind=i*nz+k;
         x_ind=i*ny*nz+(ny-1)*nz+k;
         sendBC[ind]=data[x_ind];
         x_ind=i*ny*nz+k;
         sendAD[ind]=data[x_ind];
     //    ind++;
     }
    }
    
    ind=0;  x_ind=0;
    for(k=0;k<nz;k++)
    { 
      x_ind=k;
      sendA[ind]=data[x_ind];
      x_ind=(ny-1)*nz+k;
      sendB[ind]=data[x_ind];
      x_ind=(nx-1)*ny*nz+(ny-1)*nz+k;
      sendC[ind]=data[x_ind];
      x_ind=(nx-1)*ny*nz+k;
      sendD[ind]=data[x_ind];
      ind++;
    }
}

//communite
void CommData(double *data,Info *info)
{

    int myrank=info->myrank;
    int s_ab,s_bc,s_cd ,s_ad,s_a,s_b,s_c,s_d;
    int nx,ny,nz,n;
  
    nx=info->nx;        ny=info->ny;    nz=info->nz;
    s_ab=N;   s_bc=s_ab+ny*nz;
    s_cd=s_bc+nx*nz;
    s_ad=s_cd+ny*nz;
    s_a=s_ad+nx*nz;
    s_b=s_a+nz;
    s_c=s_b+nz;
    s_d=s_c+nz;

    double *sendAB,*sendBC,*sendCD,*sendAD,*sendA,*sendB,*sendC,*sendD;
    int abSize=ny*nz;
    int bcSize=nx*nz;
    int zSize=nz;

    sendAB=(double*)malloc(abSize*sizeof(double));  sendA=(double*)malloc(zSize*sizeof(double));
    sendBC=(double*)malloc(bcSize*sizeof(double));  sendB=(double*)malloc(zSize*sizeof(double));
    sendCD=(double*)malloc(abSize*sizeof(double));  sendC=(double*)malloc(zSize*sizeof(double));
    sendAD=(double*)malloc(bcSize*sizeof(double));  sendD=(double*)malloc(zSize*sizeof(double));
    

     int reqCount=16;
     MPI_Status m_status[16];  
     MPI_Request req[16]; 
     GetSendData(sendAB,sendBC,sendCD,sendAD,sendA,sendB,sendC,sendD,data,info);  
      switch(myrank){
          case 0:
          {
            MPI_Isend(&sendAB[0],abSize,MPI_DOUBLE,6,0,MPI_COMM_WORLD,&req[0]);
            MPI_Isend(&sendBC[0],bcSize,MPI_DOUBLE,1,0,MPI_COMM_WORLD,&req[1]);
            MPI_Isend(&sendCD[0],abSize,MPI_DOUBLE,2,0,MPI_COMM_WORLD,&req[2]);
            MPI_Isend(&sendAD[0],bcSize,MPI_DOUBLE,4,0,MPI_COMM_WORLD,&req[3]);

            MPI_Isend(&sendA[0],zSize,MPI_DOUBLE,2,10,MPI_COMM_WORLD,&req[4]);
            MPI_Isend(&sendB[0],zSize,MPI_DOUBLE,7,10,MPI_COMM_WORLD,&req[5]);
            MPI_Isend(&sendC[0],zSize,MPI_DOUBLE,3,10,MPI_COMM_WORLD,&req[6]);
            MPI_Isend(&sendD[0],zSize,MPI_DOUBLE,6,10,MPI_COMM_WORLD,&req[7]);
           
      //rev 
      MPI_Irecv(&data[s_ab],abSize,MPI_DOUBLE,6,6,MPI_COMM_WORLD,&req[8]);
      MPI_Irecv(&data[s_bc],bcSize,MPI_DOUBLE,1,1,MPI_COMM_WORLD,&req[9]);
      MPI_Irecv(&data[s_cd],abSize,MPI_DOUBLE,2,2,MPI_COMM_WORLD,&req[10]);
            MPI_Irecv(&data[s_ad],bcSize,MPI_DOUBLE,4,4,MPI_COMM_WORLD,&req[11]);
            
      MPI_Irecv(&data[s_a],zSize,MPI_DOUBLE,2,12,MPI_COMM_WORLD,&req[12]);
            MPI_Irecv(&data[s_b],zSize,MPI_DOUBLE,7,17,MPI_COMM_WORLD,&req[13]);
            MPI_Irecv(&data[s_c],zSize,MPI_DOUBLE,3,13,MPI_COMM_WORLD,&req[14]);
            MPI_Irecv(&data[s_d],zSize,MPI_DOUBLE,6,16,MPI_COMM_WORLD,&req[15]);
            break;
          }
          case 1:
          {
            MPI_Isend(&sendAB[0],abSize,MPI_DOUBLE,7,1,MPI_COMM_WORLD,&req[0]);
            MPI_Isend(&sendBC[0],bcSize,MPI_DOUBLE,5,1,MPI_COMM_WORLD,&req[1]);
            MPI_Isend(&sendCD[0],abSize,MPI_DOUBLE,3,1,MPI_COMM_WORLD,&req[2]);
            MPI_Isend(&sendAD[0],bcSize,MPI_DOUBLE,0,1,MPI_COMM_WORLD,&req[3]);

            MPI_Isend(&sendA[0],zSize,MPI_DOUBLE,6,11,MPI_COMM_WORLD,&req[4]);
            MPI_Isend(&sendB[0],zSize,MPI_DOUBLE,3,11,MPI_COMM_WORLD,&req[5]);
            MPI_Isend(&sendC[0],zSize,MPI_DOUBLE,7,11,MPI_COMM_WORLD,&req[6]);
            MPI_Isend(&sendD[0],zSize,MPI_DOUBLE,2,11,MPI_COMM_WORLD,&req[7]);
           
      //rev 
      MPI_Irecv(&data[s_ab],abSize,MPI_DOUBLE,7,7,MPI_COMM_WORLD,&req[8]);
      MPI_Irecv(&data[s_bc],bcSize,MPI_DOUBLE,5,5,MPI_COMM_WORLD,&req[9]);
      MPI_Irecv(&data[s_cd],abSize,MPI_DOUBLE,3,3,MPI_COMM_WORLD,&req[10]);
            MPI_Irecv(&data[s_ad],bcSize,MPI_DOUBLE,0,0,MPI_COMM_WORLD,&req[11]);
            
      MPI_Irecv(&data[s_a],zSize,MPI_DOUBLE,6,16,MPI_COMM_WORLD,&req[12]);
            MPI_Irecv(&data[s_b],zSize,MPI_DOUBLE,3,13,MPI_COMM_WORLD,&req[13]);
            MPI_Irecv(&data[s_c],zSize,MPI_DOUBLE,7,17,MPI_COMM_WORLD,&req[14]);
            MPI_Irecv(&data[s_d],zSize,MPI_DOUBLE,2,12,MPI_COMM_WORLD,&req[15]);
            break;
          }
          case 2:
          {
            MPI_Isend(&sendAB[0],abSize,MPI_DOUBLE,0,2,MPI_COMM_WORLD,&req[0]);
            MPI_Isend(&sendBC[0],bcSize,MPI_DOUBLE,3,2,MPI_COMM_WORLD,&req[1]);
            MPI_Isend(&sendCD[0],abSize,MPI_DOUBLE,4,2,MPI_COMM_WORLD,&req[2]);
            MPI_Isend(&sendAD[0],bcSize,MPI_DOUBLE,6,2,MPI_COMM_WORLD,&req[3]);

            MPI_Isend(&sendA[0],zSize,MPI_DOUBLE,4,12,MPI_COMM_WORLD,&req[4]);
            MPI_Isend(&sendB[0],zSize,MPI_DOUBLE,1,12,MPI_COMM_WORLD,&req[5]);
            MPI_Isend(&sendC[0],zSize,MPI_DOUBLE,5,12,MPI_COMM_WORLD,&req[6]);
            MPI_Isend(&sendD[0],zSize,MPI_DOUBLE,0,12,MPI_COMM_WORLD,&req[7]);
           
      //rev 
      MPI_Irecv(&data[s_ab],abSize,MPI_DOUBLE,0,0,MPI_COMM_WORLD,&req[8]);
      MPI_Irecv(&data[s_bc],bcSize,MPI_DOUBLE,3,3,MPI_COMM_WORLD,&req[9]);
      MPI_Irecv(&data[s_cd],abSize,MPI_DOUBLE,4,4,MPI_COMM_WORLD,&req[10]);
            MPI_Irecv(&data[s_ad],bcSize,MPI_DOUBLE,6,6,MPI_COMM_WORLD,&req[11]);
            
      MPI_Irecv(&data[s_a],zSize,MPI_DOUBLE,4,14,MPI_COMM_WORLD,&req[12]);
            MPI_Irecv(&data[s_b],zSize,MPI_DOUBLE,1,11,MPI_COMM_WORLD,&req[13]);
            MPI_Irecv(&data[s_c],zSize,MPI_DOUBLE,5,15,MPI_COMM_WORLD,&req[14]);
            MPI_Irecv(&data[s_d],zSize,MPI_DOUBLE,0,10,MPI_COMM_WORLD,&req[15]);
            break;
          }
          case 3:
          {
            MPI_Isend(&sendAB[0],abSize,MPI_DOUBLE,1,3,MPI_COMM_WORLD,&req[0]);
            MPI_Isend(&sendBC[0],bcSize,MPI_DOUBLE,7,3,MPI_COMM_WORLD,&req[1]);
            MPI_Isend(&sendCD[0],abSize,MPI_DOUBLE,5,3,MPI_COMM_WORLD,&req[2]);
            MPI_Isend(&sendAD[0],bcSize,MPI_DOUBLE,2,3,MPI_COMM_WORLD,&req[3]);

            MPI_Isend(&sendA[0],zSize,MPI_DOUBLE,0,13,MPI_COMM_WORLD,&req[4]);
            MPI_Isend(&sendB[0],zSize,MPI_DOUBLE,5,13,MPI_COMM_WORLD,&req[5]);
            MPI_Isend(&sendC[0],zSize,MPI_DOUBLE,1,13,MPI_COMM_WORLD,&req[6]);
            MPI_Isend(&sendD[0],zSize,MPI_DOUBLE,4,13,MPI_COMM_WORLD,&req[7]);
           
      //rev 
      MPI_Irecv(&data[s_ab],abSize,MPI_DOUBLE,1,1,MPI_COMM_WORLD,&req[8]);
      MPI_Irecv(&data[s_bc],bcSize,MPI_DOUBLE,7,7,MPI_COMM_WORLD,&req[9]);
      MPI_Irecv(&data[s_cd],abSize,MPI_DOUBLE,5,5,MPI_COMM_WORLD,&req[10]);
            MPI_Irecv(&data[s_ad],bcSize,MPI_DOUBLE,2,2,MPI_COMM_WORLD,&req[11]);
            
      MPI_Irecv(&data[s_a],zSize,MPI_DOUBLE,0,10,MPI_COMM_WORLD,&req[12]);
            MPI_Irecv(&data[s_b],zSize,MPI_DOUBLE,5,15,MPI_COMM_WORLD,&req[13]);
            MPI_Irecv(&data[s_c],zSize,MPI_DOUBLE,1,11,MPI_COMM_WORLD,&req[14]);
            MPI_Irecv(&data[s_d],zSize,MPI_DOUBLE,4,14,MPI_COMM_WORLD,&req[15]);
            break;
          }
          case 4:
          {
            MPI_Isend(&sendAB[0],abSize,MPI_DOUBLE,2,4,MPI_COMM_WORLD,&req[0]);
            MPI_Isend(&sendBC[0],bcSize,MPI_DOUBLE,5,4,MPI_COMM_WORLD,&req[1]);
            MPI_Isend(&sendCD[0],abSize,MPI_DOUBLE,6,4,MPI_COMM_WORLD,&req[2]);
            MPI_Isend(&sendAD[0],bcSize,MPI_DOUBLE,0,4,MPI_COMM_WORLD,&req[3]);

            MPI_Isend(&sendA[0],zSize,MPI_DOUBLE,6,14,MPI_COMM_WORLD,&req[4]);
            MPI_Isend(&sendB[0],zSize,MPI_DOUBLE,3,14,MPI_COMM_WORLD,&req[5]);
            MPI_Isend(&sendC[0],zSize,MPI_DOUBLE,7,14,MPI_COMM_WORLD,&req[6]);
            MPI_Isend(&sendD[0],zSize,MPI_DOUBLE,2,14,MPI_COMM_WORLD,&req[7]);
           
      //rev 
      MPI_Irecv(&data[s_ab],abSize,MPI_DOUBLE,2,2,MPI_COMM_WORLD,&req[8]);
      MPI_Irecv(&data[s_bc],bcSize,MPI_DOUBLE,5,5,MPI_COMM_WORLD,&req[9]);
      MPI_Irecv(&data[s_cd],abSize,MPI_DOUBLE,6,6,MPI_COMM_WORLD,&req[10]);
            MPI_Irecv(&data[s_ad],bcSize,MPI_DOUBLE,0,0,MPI_COMM_WORLD,&req[11]);
            
      MPI_Irecv(&data[s_a],zSize,MPI_DOUBLE,6,16,MPI_COMM_WORLD,&req[12]);
            MPI_Irecv(&data[s_b],zSize,MPI_DOUBLE,3,13,MPI_COMM_WORLD,&req[13]);
            MPI_Irecv(&data[s_c],zSize,MPI_DOUBLE,7,17,MPI_COMM_WORLD,&req[14]);
            MPI_Irecv(&data[s_d],zSize,MPI_DOUBLE,2,12,MPI_COMM_WORLD,&req[15]);
            break;
          }
          case 5:
          {
            MPI_Isend(&sendAB[0],abSize,MPI_DOUBLE,3,5,MPI_COMM_WORLD,&req[0]);
            MPI_Isend(&sendBC[0],bcSize,MPI_DOUBLE,1,5,MPI_COMM_WORLD,&req[1]);
            MPI_Isend(&sendCD[0],abSize,MPI_DOUBLE,7,5,MPI_COMM_WORLD,&req[2]);
            MPI_Isend(&sendAD[0],bcSize,MPI_DOUBLE,4,5,MPI_COMM_WORLD,&req[3]);

            MPI_Isend(&sendA[0],zSize,MPI_DOUBLE,2,15,MPI_COMM_WORLD,&req[4]);
            MPI_Isend(&sendB[0],zSize,MPI_DOUBLE,7,15,MPI_COMM_WORLD,&req[5]);
            MPI_Isend(&sendC[0],zSize,MPI_DOUBLE,3,15,MPI_COMM_WORLD,&req[6]);
            MPI_Isend(&sendD[0],zSize,MPI_DOUBLE,6,15,MPI_COMM_WORLD,&req[7]);
           
      //rev 
      MPI_Irecv(&data[s_ab],abSize,MPI_DOUBLE,3,3,MPI_COMM_WORLD,&req[8]);
      MPI_Irecv(&data[s_bc],bcSize,MPI_DOUBLE,1,1,MPI_COMM_WORLD,&req[9]);
      MPI_Irecv(&data[s_cd],abSize,MPI_DOUBLE,7,7,MPI_COMM_WORLD,&req[10]);
            MPI_Irecv(&data[s_ad],bcSize,MPI_DOUBLE,4,4,MPI_COMM_WORLD,&req[11]);
            
      MPI_Irecv(&data[s_a],zSize,MPI_DOUBLE,2,12,MPI_COMM_WORLD,&req[12]);
            MPI_Irecv(&data[s_b],zSize,MPI_DOUBLE,7,17,MPI_COMM_WORLD,&req[13]);
            MPI_Irecv(&data[s_c],zSize,MPI_DOUBLE,3,13,MPI_COMM_WORLD,&req[14]);
            MPI_Irecv(&data[s_d],zSize,MPI_DOUBLE,6,16,MPI_COMM_WORLD,&req[15]);
            break;
          }
          case 6:
          {
            MPI_Isend(&sendAB[0],abSize,MPI_DOUBLE,4,6,MPI_COMM_WORLD,&req[0]);
            MPI_Isend(&sendBC[0],bcSize,MPI_DOUBLE,7,6,MPI_COMM_WORLD,&req[1]);
            MPI_Isend(&sendCD[0],abSize,MPI_DOUBLE,0,6,MPI_COMM_WORLD,&req[2]);
            MPI_Isend(&sendAD[0],bcSize,MPI_DOUBLE,2,6,MPI_COMM_WORLD,&req[3]);

            MPI_Isend(&sendA[0],zSize,MPI_DOUBLE,0,16,MPI_COMM_WORLD,&req[4]);
            MPI_Isend(&sendB[0],zSize,MPI_DOUBLE,5,16,MPI_COMM_WORLD,&req[5]);
            MPI_Isend(&sendC[0],zSize,MPI_DOUBLE,1,16,MPI_COMM_WORLD,&req[6]);
            MPI_Isend(&sendD[0],zSize,MPI_DOUBLE,4,16,MPI_COMM_WORLD,&req[7]);
           
      //rev 
      MPI_Irecv(&data[s_ab],abSize,MPI_DOUBLE,4,4,MPI_COMM_WORLD,&req[8]);
      MPI_Irecv(&data[s_bc],bcSize,MPI_DOUBLE,7,7,MPI_COMM_WORLD,&req[9]);
      MPI_Irecv(&data[s_cd],abSize,MPI_DOUBLE,0,0,MPI_COMM_WORLD,&req[10]);
            MPI_Irecv(&data[s_ad],bcSize,MPI_DOUBLE,2,2,MPI_COMM_WORLD,&req[11]);
            
      MPI_Irecv(&data[s_a],zSize,MPI_DOUBLE,0,10,MPI_COMM_WORLD,&req[12]);
            MPI_Irecv(&data[s_b],zSize,MPI_DOUBLE,5,15,MPI_COMM_WORLD,&req[13]);
            MPI_Irecv(&data[s_c],zSize,MPI_DOUBLE,1,11,MPI_COMM_WORLD,&req[14]);
            MPI_Irecv(&data[s_d],zSize,MPI_DOUBLE,4,14,MPI_COMM_WORLD,&req[15]);
            break;
          }
          case 7:
          {
            MPI_Isend(&sendAB[0],abSize,MPI_DOUBLE,5,7,MPI_COMM_WORLD,&req[0]);
            MPI_Isend(&sendBC[0],bcSize,MPI_DOUBLE,3,7,MPI_COMM_WORLD,&req[1]);
            MPI_Isend(&sendCD[0],abSize,MPI_DOUBLE,1,7,MPI_COMM_WORLD,&req[2]);
            MPI_Isend(&sendAD[0],bcSize,MPI_DOUBLE,6,7,MPI_COMM_WORLD,&req[3]);

            MPI_Isend(&sendA[0],zSize,MPI_DOUBLE,4,17,MPI_COMM_WORLD,&req[4]);
            MPI_Isend(&sendB[0],zSize,MPI_DOUBLE,1,17,MPI_COMM_WORLD,&req[5]);
            MPI_Isend(&sendC[0],zSize,MPI_DOUBLE,5,17,MPI_COMM_WORLD,&req[6]);
            MPI_Isend(&sendD[0],zSize,MPI_DOUBLE,0,17,MPI_COMM_WORLD,&req[7]);
           
      //rev 
      MPI_Irecv(&data[s_ab],abSize,MPI_DOUBLE,5,5,MPI_COMM_WORLD,&req[8]);
      MPI_Irecv(&data[s_bc],bcSize,MPI_DOUBLE,3,3,MPI_COMM_WORLD,&req[9]);
      MPI_Irecv(&data[s_cd],abSize,MPI_DOUBLE,1,1,MPI_COMM_WORLD,&req[10]);
            MPI_Irecv(&data[s_ad],bcSize,MPI_DOUBLE,6,6,MPI_COMM_WORLD,&req[11]);
            
      MPI_Irecv(&data[s_a],zSize,MPI_DOUBLE,4,14,MPI_COMM_WORLD,&req[12]);
            MPI_Irecv(&data[s_b],zSize,MPI_DOUBLE,1,11,MPI_COMM_WORLD,&req[13]);
            MPI_Irecv(&data[s_c],zSize,MPI_DOUBLE,5,15,MPI_COMM_WORLD,&req[14]);
            MPI_Irecv(&data[s_d],zSize,MPI_DOUBLE,0,10,MPI_COMM_WORLD,&req[15]);
            break;
          }
          default:
       break;
  }
      MPI_Waitall(reqCount, req, m_status);  
      free(sendAB);free(sendBC); free(sendCD);  free(sendAD);  free(sendA);  free(sendB); free(sendC); free(sendD); 
}

//communication
void CGRCommunite(Info *info)
{
    int myrank=info->myrank;
    int ind=0;
    int i,j;
    double *R=(double*)malloc(nN*sizeof(double));
    double *Rm=(double*)malloc(nN*sizeof(double));
    double *p=(double*)malloc(nN*sizeof(double));
    double *Ap=(double*)malloc(nN*sizeof(double));
    double *ARm=(double*)malloc(nN*sizeof(double));
   // double *tmp=(double*)malloc(nN*sizeof(double));

    double *Apj[5];
    double *pj[5];
    double betaij[5];
    double deta;
    int k,m_ind;
    for(i=0;i<5;i++){
        pj[i]=(double*)malloc(nN*sizeof(double));
        Apj[i]=(double*)malloc(nN*sizeof(double));
        betaij[i]=0;
    }
    ind=0;
    getAmVector(R,x0);
    #pragma omp parallel for private(i) shared(R,b)  
    for(i=0;i<N;i++){
       R[i]=b[i]-R[i];
    }
 
    CommData(R,info);
    JacobiCal(Rm,R); 
   //  getMatrixLU(Rm,R); 
     memcpy(p,Rm,N*sizeof(double));

    double ai,checkVal;
    double revSum;
    double sendVal; 
    int m=0;
    i=1; k=5;    //子空间
    do{
  CommData(p,info);
  
        getAmVector(Ap,p);   //get Api-1 向量 
    memcpy(pj[(i-1)%k],p,N*sizeof(double));
    memcpy(Apj[(i-1)%k],Ap,N*sizeof(double));

        sendVal=VectorDianJi(Ap,Ap);
        MPI_Allreduce(&sendVal,&revSum,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
        betaij[(i-1)%k]=revSum;      

        sendVal=VectorDianJi(R,Ap);
        MPI_Allreduce(&sendVal,&revSum,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
  ai=revSum/betaij[(i-1)%k];

  #pragma omp parallel for private(j) shared(x0,R,Ap,ai,p)  
        for(j=0;j<N;j++){       //xi=x(i-1)+a(i-1)*p(i-1)
            x0[j]=x0[j]+ai*p[j];
            R[j]=R[j]-ai*Ap[j];
//      Rm[j]=R[j];  
        }

        CommData(R,info);
        JacobiCal(Rm,R);        //预处理
       // getMatrixLU(Rm,R);       

  sendVal=VectorDianJi(Rm,Rm);
  MPI_Allreduce(&sendVal,&revSum,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
//  if(myrank==0)printf("checkSum---->%.30lf \n",revSum); 
  if(revSum<0.0000000001) break;

       CommData(Rm,info);
       memcpy(p,Rm,N*sizeof(double));

       getAmVector(ARm,Rm);
       ind=0;
       for(j=((i-1)/k)*k;j<i;j++){
     sendVal=VectorDianJi(ARm,Apj[ind]);
           MPI_Allreduce(&sendVal,&revSum,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); 
           deta=-revSum/betaij[ind];
     #pragma omp parallel for private(m) shared(deta,pj)   
           for(m=0;m<N;m++){
             p[m]=p[m]+deta*pj[ind][m];
           }
    ind++;   
        }       
        i++;  
     }while(i<150);
     if(myrank==0)   printf("Jacobi times------>%d\n",(i-1));
 
    for(i=0;i<5;i++){
        free(pj[i]);  free(Apj[i]);
     }
     free(R);free(Rm);free(p);free(Ap);free(ARm);   
}

double difftimeval(const struct timeval *start, const struct timeval *end)
{
        double d;
        time_t s;
        suseconds_t u;
        s = start->tv_sec - end->tv_sec;
        u = start->tv_usec - end->tv_usec;
        d = s;
        d += u/1000000.0;
       return d;
}

int main(int argc, char **argv) {
   AJ=NULL; AV=NULL; b=NULL; x0=NULL;
   MPI_Init(&argc,&argv);
   int myrank;
   int NX,NY,NZ;
   int PX,PY,PZ;
   NX=atoi(argv[1]);NY = atoi(argv[2]); NZ = atoi(argv[3]); 
   PX=atoi(argv[4]);PY = atoi(argv[5]); PZ = atoi(argv[6]); 
   char *filename_A=argv[7];
   char *filename_x0=argv[8];
   char *filename_b=argv[9];
   if(PZ!=1) return 1;
   struct timeval t1, t2;   
   MPI_Comm_rank(MPI_COMM_WORLD,&myrank);

    nSize=NX*NY*NZ;
    Info info=InitInfo(NX,NY,NZ,PX,PY,PZ);
    MPI_Barrier(MPI_COMM_WORLD);
    gettimeofday(&t1, NULL);
    LoadData(&info,filename_A,filename_x0,filename_b);
  //  getLUmatrixA();
  //
    CGRCommunite(&info); 

    MPI_Barrier(MPI_COMM_WORLD);
    gettimeofday(&t2, NULL);
 
    if(myrank==0)printf("Total time: %.6lf\n",difftimeval(&t2,&t1));   
    free(AJ);  free(AV);  free(b);  free(x0);  free(diag);
    MPI_Finalize();
    return 0;
}






























